Instructions to access sf object mapping from Chris’ GitHub: https://github.com/cfree14/wcfish

# Run if you don't already have devtools installed
# install.packages("devtools")

# Run once devtools is successfully installed
# devtools::install_github("cfree14/wcfish", force=T)
# library(wcfish)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ tibble  3.1.8     ✔ dplyr   1.1.0
## ✔ tidyr   1.2.1     ✔ stringr 1.4.1
## ✔ readr   2.1.3     ✔ forcats 0.5.2
## ✔ purrr   0.3.5     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks plotly::filter(), stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(here)
## here() starts at /Users/clairegonzales/Documents/R projects/chap_2_sitingstudy/chap2
library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
#library(tmap)
# library(paletter)

Getting the blocks shapefile from Chris’ Github:

# See dataset metadata by typing a question mark and the dataset name
 # ?wcfish::blocks

# Store the dataset in an object of your own choosing as follows:
blocks_sf <- wcfish::blocks %>%
  drop_na(block_lat_dd, block_long_dd) 
# %>%
#   st_as_sf(coords = c("block_long_dd", "block_lat_dd")) #take lat and long and convert to spatial points

################################################################################
# FOR ACTUAL MAPPING
# Use the ESRI .shp file from: https://chrismfree.com/california-fisheries-data/ca-fishing-blocks/

blocks_shp <- read_sf(here("data_fromChris", "CA_OR_WA_comm_fishing_blocks", "CA_OR_WA_comm_fishing_blocks.shp"))

Read in catch data from Chris

catchdata <- read_csv(here("data_fromChris", "gonzalez_data.csv"))
## New names:
## Rows: 640 Columns: 8
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (1): confidential dbl (7): ...1, block_id, nvessels, nfishers, nbusinesses,
## value_usd, landing...
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`

Now that we have both the catch data from Chris and the blocks shapefile, I am going to run these through a spatial function

# looking at blocks spatially
ggplot(data = blocks_sf, aes(x = block_long_dd, y = block_lat_dd)) +
  geom_point(shape = 15, fill = NA, aes(size = block_sqkm)) +
  theme_minimal()

Making the catch data ready for spatial mapping

catchdata_sf_merge <- merge(blocks_sf, catchdata) %>% 
  filter(block_state == "California",
         block_type %in% c("Inshore", "Midshore")) %>% 
  mutate() %>% 
  select(block_id, block_type, block_sqkm, block_long_dd, block_lat_dd, nvessels:geometry) %>% 
  drop_na(block_long_dd, block_lat_dd) %>% 
  st_as_sf(coords = c("block_long_dd", "block_lat_dd")) #take lat and long and convert to spatial points

Coming back later to add a column of the z transformed data in R. This seems easier to do here than in ArcPro.

a <- min(catchdata_sf_merge$landings_lb, na.rm = TRUE)
b <- max(catchdata_sf_merge$landings_lb, na.rm = TRUE) + 1
m <- (a + b)/2
# median(catchdata_sf_merge$landings_lb, na.rm = TRUE)

catchdata_sf_merge_z <- catchdata_sf_merge %>% 
  mutate(z_lndngs_ = case_when(
    landings_lb <= a ~ 1,
    (landings_lb > a & landings_lb <= m) ~ (1-(2*(((landings_lb-a)/(b-a))**2))),
    (landings_lb > m & landings_lb < b) ~ (2*(((landings_lb - b)/(b-a))**2)),
    landings_lb >= b ~ 0,
    TRUE ~ as.numeric(NA)
  ))

# trying to normalize the data by taking the natural log
catchdata_sf_log <- catchdata_sf_merge %>% 
  mutate(ln_landngs = log(landings_lb))
# WORKED. But also now everything is negative soooo let's see how that looks

Trying a linear transformation of the data

slope <- (1/-(b-a))
slope2 <- (1/-b)
y_int <- 1

catchdata_sf_merge_linear <- catchdata_sf_merge %>% 
  mutate(linear_landings = (landings_lb*slope2)+1)

## See bottom chunk for plots and plotly of this transformed data

Option 4 - log transformation of data and then Z transformation

a_ln <- min(catchdata_sf_log$ln_landngs, na.rm = TRUE)
b_ln <- max(catchdata_sf_log$ln_landngs, na.rm = TRUE) + 1
m_ln <- (a_ln + b_ln)/2

catchdata_sf_merge_ln_z <- catchdata_sf_log %>% 
  mutate(z_lndngs_ln = case_when(
    ln_landngs <= a_ln ~ 1,
    (ln_landngs > a_ln & ln_landngs <= m_ln) ~ (1-(2*(((ln_landngs-a_ln)/(b_ln-a_ln))**2))),
    (ln_landngs > m_ln & ln_landngs < b_ln) ~ (2*(((ln_landngs - b_ln)/(b_ln-a_ln))**2)),
    ln_landngs >= b_ln ~ 0,
    TRUE ~ as.numeric(NA)
  ))

Let’s take a look at the transformed data

#scatter plot of ln and z transformed data, against original catch data
ggplot(data = catchdata_sf_log, aes(x = landings_lb, y = ln_landngs)) +
  geom_point() +
  theme_minimal()
## Warning: Removed 1 rows containing missing values (`geom_point()`).

# histogram of original catch data
ggplot(data = catchdata_sf_log, aes(x = landings_lb)) +
  geom_histogram() +
  theme_minimal() +
  xlim(0, 2500000)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 27 rows containing non-finite values (`stat_bin()`).
## Warning: Removed 2 rows containing missing values (`geom_bar()`).

# histogram of ln and z transformed data
ggplot(data = catchdata_sf_log, aes(x = ln_landngs)) +
  geom_histogram() +
  theme_minimal()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

# we need to tell r what the coord reference system is using `st_crs`

#st_crs(catchdata_sf_merge) <- 4326

#UPDATE
# Doing this with the transformed data instead

st_crs(catchdata_sf_merge_z) <- 4326

Add blocks shapefile to map

# 
# blocks_shp_transform <- st_transform(blocks_shp, 4326)
# 
# ggplot(data = blocks_shp_transform) +
#   geom_sf()+
#   theme_minimal()

Combine catchdata and block shapfile map

# map <- ggplot() +
#   geom_sf(data = blocks_shp_transform) +
#   geom_sf(data = catchdata_sf_merge, na.rm = TRUE,
#           aes(color = landings_lb, size = block_sqkm)) +
#   theme_void()

#trying this a different way

catchdata_sf_merge_transform <- st_transform(catchdata_sf_merge_z, 4326)

#map of landings
ggplot(data = catchdata_sf_merge_transform) +
  geom_sf(aes(fill = landings_lb)) +
  guides(color = guide_legend(reverse = TRUE)) +
  theme_minimal()

#map of value
ggplot(data = catchdata_sf_merge_transform) +
  geom_sf(aes(fill = value_usd))+
  theme_minimal()

# looking at data

ggplot(data = catchdata_sf_merge_transform, aes(y = landings_lb)) +
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

# Looks like the data is strongly skewed left, which is why the legend is so not useful

Write the catch data as a shapefile so that we can upload to ArcPro

st_write(catchdata_sf_merge_transform, here("data", "catchdata_output", "catchdata_output.shp"), append = FALSE)
## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver
## Deleting layer `catchdata_output' using driver `ESRI Shapefile'
## Writing layer `catchdata_output' to data source 
##   `/Users/clairegonzales/Documents/R projects/chap_2_sitingstudy/chap2/data/catchdata_output/catchdata_output.shp' using driver `ESRI Shapefile'
## Writing 547 features with 12 fields and geometry type Multi Polygon.

Visualizations

OPTIONS: Data that is transformed to Z-shaped membership function OR Linear transformation

plotz <- ggplot(data = catchdata_sf_merge_z, aes(x = landings_lb, y = z_lndngs_)) +
  geom_point() +
  theme_minimal()

plotz_plotly <- ggplotly(plotz)

plotz_plotly
#scatter plot of ln and z transformed data, against original catch data
plot1 <- ggplot(data = catchdata_sf_log, aes(x = landings_lb, y = ln_landngs)) +
  geom_point() +
  theme_minimal()

plot1_plotly <- ggplotly(plot1)

plot1_plotly
plot2 <- ggplot(data = catchdata_sf_merge_linear, aes(x = landings_lb, y = linear_landings)) +
  geom_point() +
  theme_minimal()

plot2_plotly <- ggplotly(plot2)

plot2_plotly
plot4 <- ggplot(data = catchdata_sf_merge_ln_z, aes(x = ln_landngs, y = z_lndngs_ln)) +
  geom_point() +
  theme_minimal()

plot4_plotly <- ggplotly(plot4)

plot4_plotly